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© Digital image processing taking into account the standard deviation. 

© Image processing methods are described that 
take into account digitized pixel values, the average 
pixel value within a window surrounding a pixel of 
interest and the standard deviation, jol-.sajd, pjxel 
values„and that can be performed with less com- 
putational effort than the prior art methods. 
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Digital image processing taking into account the standard deviation 



BACKGROUND OF THE INVENTION 



Field of the invention 

The present invention is in the field of digital 
image processing. The invention more particulary 
relates to digital image processing to be used for 
the purpose of image enhancement in the field of 
computer assisted diagnostic imaging systems. 



Description of the prior art 

A general problem in radiology is the discer- 
nability of small details especially if their density is 
quite similar to that of the surrounding area. 
Image processing methods were developed by 
means of which the detectability of details can be 
enhanced and by means of which the information 
transfer from an image to the observer can be 
maximized. By means of said image processing 
methods e.g. edge .enhancement and contrast en- 
- hancement can be obtained. - 
image^rolTe^"ncpm~radiology can be performed 
on a digital representation of a radiologic image 
that is obtained either through a radiologic method 
in which a digital image of a body is directly 
obtained, such as computerized tomography, or 
through a radiologic method in which a digital im- 
age is indirectly obtained e.g. by scanning an im- 
age stored in an intermediate storage device such 
as a silver halide emulsion film or a photostimula- 
ble phosphor. 

In prior art image processing methods, e.g. the 
unsharp masking method, the resulting image after 
processing is function of the average pixel value. 
This pixel value can be a density value, an inten- 
sity value or the logarithm of an intensity value. In 
the following whenever a pixel density value is 
used this could be replaced by an intensity value 
or by the logarithm of an intensity value. 
However, the above-mentioned average pixel value 
gives poor indication of the local image information. 
A better indicator is the standard deviation of the 
pixel densities which factor gives an indication of 
the pixel density distribution in a region surround- 
ing the pixel of interest.By using the standard de- 
viation, the image enhancement factor can be con- 
trolled locally so that the noise in an image can be 
taken into account. 

Prior art image processing methods that take 
into account the standard deviation of the pixels as 
well as the average pixel value and the original 
pixel values are known and are e.g. described in a 



publication by M.I. Cocklin et al. in Image and 
Vision Computing, Vol. 1, Nr. 2. May 1983, p. 70 
and following pages. 

Although a lot of research has been performed in 
5 respect of image processing methods that take into 
account the standard deviation, a prejudice still 
exists against the use of the standard deviation 
because of the increase of computational effort that 
this usually entails.Therefore several prior art ap- 
70 proaches use an approximation of the standard 
deviation of the pixel values. 

Among existing image processing techniques 
adaptive histogram equalisation is in respect of the 
present invention of particular interest. 

1 75 Histogram equalisation is an image processing 
technique by means of which overall contrast is 
improved by equally distributing the pixel densities 
in an image throughout the available grey scale. 
According to this technique the histogram and a 
20 rescaled cumulative histogram of an image are 
calculated and each pixel of an original image is 
mapped to its corresponding value in said rescaled 

cumulative histogram. . 

For some pictures better results (improved contrast 

25 in smaller regions) are obtained by applying the 
continuous adaptive histogram equalisation tech- 
nique according to which an image is scanned by 
means of a sliding window (a region of N times M 
pixels surrounding a pixel of interest and sliding 

30 over the entire image). Again the histogram and a 
rescaled cumulative histogram is calculated for 
each area covered by the sliding window. The 
central pixel in the window is mapped to its cor- 
responding value in the rescaled cumulative histo- 

35 gram. This procedure is repeated for each pixel. 

A variant is the adaptive histogram equalisation 
technique according to which an image is sub- 
divided into regions or windows. 
The histogram and the cumulative histogram are 

40 calculated for each window and the pixel of interest 
is mapped to its corresponding value in said cu- 
mulative histogram. Smooth transitions between in- 
dividual zones can for example be obtained by 
means of interpolation techniques. 

45 To prevent the overamplification in regions with a 
very narrow histogram a limit can be set on the 
maximum slope of the cumulative histogram for 
example by clipping the original histogram of the 
image. This method is commonly called contrast 

so limited adaptive histogram equalisation. 

An overview of the above techniques is i.e. given in 
the introductory part of a publication by R.H. Sher- 
rier and G.A. Johnson in IEEE Transactions on 
Medical Imaging Vol. M1-6, No. 1, March 1987. 
The above methods have some disadvantages: the 
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methods are rather slow and ask for a lot of com- 
putational effort and under certain conditions the 
processed image shows rather undesirable effects. 
An amended form, of the histogram equalisation 
technique in which the standard deviation of the 
histogram distribution of a region surrounding a 
pixel cf interest is taken into account is disclosed in 
a publication by Hartwig Blume and Kanji Kamiya 
in SPIE Vol. 767 Medical Imaging (1987) p. 375. 
In the captioned disclosure is published that such 
procedure would greatly increase the computation- 
al effort. Hence a prejudice exists against the use 
of the local standard deviation in an amended form 
of adaptive histogram equalisation technique. 
A method for exerting the above form of histogram 
equalisation in which the standard deviation is tak- 
en into account within a reasonable amount of 
computation time is not known . 

Another existing image processing technique 
namely the unsharp masking technique is in re- 
spect of the present invention of particular inter- 
est.Commonly the result of the image processing 
technique can be formulated as follows: y = x + 
alfa.T(x-m) wherein y represents a pixel value -after 
image processing, x represents an original pixel 
value, alfa determines the amplification, m is an 
average pixel value and T is an S-like function with 
in the normalised form the following properties T- 
(0) = 0 and t'(0) = 1 (T representing the first deri- 
vate of the function T). This function is used to 
prevent overamplification of large values x- 
m.Commonly the amplification factor is made func- 
tion of x and/or m so that an estimation is made of 
the pixel vaiues in an image region surrounding a 
pixel of interest 

The unsharp masking formula can be reformulated 
more generally as follows: y = x + alfa> (x).alfa 2 - 
(m).T(x-m). 

The pixel value nor the average pixel value 
gives a satisfactory indication of the image informa- 
tion in a region surrounding a pixel of interest since 
these values depend on the irradiation dose (x-ray 
dose) that is used. 

The standard deviation gives a better indication of 
the local image information. A variant of the un- 
sharp masking method that takes into account the 
standard deviation can be formulated for example 
as follows: y = x + alfa(s).T(x-m) wherein s is the 
stand^r^^j/iatir^n . 

One disadvantage of the unsharp masking method 
of the state of the art is the need of a smoothing 
step after application of the unsharp masking meth- 
od due to excessive noise enhancement when us- 
ing a rather large amplification factor. 
Another disadvantage is the large amount of com- 
putational effort that is needed to perform unsharp 
masking taking into account the standard deviation. 
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SUMMARY OF THE INVENTION 



Objects of the invention 

An object of the present invention is to provide 
image processing methods that take into account 
the standard deviation and that demand less com- 
putational efforts than the prior art methods. 

A further object of the present invention is to 
provide an amended form of adaptive histogram 
equalisation and contrast limited ad aptive histo- 
gra m eq ualisation- and an amended form of the 
unsharp masking that have the above enumerated 
characteristics, namely that take into account the 
standard deviation and need less computational 
effort than the prior art methods. 

Another object of the present invention is to 
provide image processing methods that take into 
account the standard deviation wherein an amend- 
ed form of standard deviation is used which takes 
less computational efforts. 

Other objects will become apparant from the 
description hereinafter. 

Image processing that takes into account the 
standard deviation is called second moment image 
processing (SMIP). 



30 Statement of invention 



According to the present invention there is 
provided an image processing method wherein a 
representation of an image in the form of digital 
pixel values x is processed to obtain resulting 
digital pixel values y as a function of a pixel value 
x, an average value m^oLpixel-vaiues-in-a-windo.w_ 
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su rrounding said pixej_x_and a standard deviation s 
of pixel vaiues in such window^ characterised "Tn 
thaf^safd' standard deviation s is obtained as ( 
S 2 /(N 2 .M 2 ) -(Si/fNLMi)) 2 ) 1/2 wherein Si is ob- 
tained as the result of a convolution of pixel values 
x and a kernel of Ni times Mi identical elements 
surrounding said pixel x and S2 is obtained as the 
result of a convolution of the squares of said pixel 
values x.said square values being obtained through 
the use of a look-up table, and a kernel of N 2 times 
M 2 identical elements surrounding said pixel x and 
wherein the square root values are obtained 
through the use of a look-up table. 

Preferably the values N1.M1 and/or N 2 .M 2 are 
a power of 2. 

The average value m in the above method is 
preferably obtained as Si/(Ni.M-,) wherein Si ,Ni 
,and Mi are as defined hereinbefore. 

In a particular application this method is used 
in diagnostic imaging methods wherein said pixels 
x are digitized density or intensity values (or the 



3 



SDOCID: <EP 0357B42A1_I_> 



5 EP 0 357 842 A1 



6 



logarithm of such intensity values) obtained by 
scanning and digitizing of a radiologic image stored 
on photographic film material or a photostimuiable 
phoshor. 

Alternatively this method can be applied in 
radiologic methods wherein no intermediary stor- 
age medium is used e.g. computer tomography.lt 
will be clear that other applications within the field 
of diagnostic imaging as well as image processing 
applications other than those related to diagnostic 
imaging may be envisioned. 

According to statistical theory the standard de- 
viation (furtheron referred to as s) of a set of N 
values denoted by "x" is obtained as (Sum (x 2 )/N - 
(Sum(x)/N) 2 ) ,/2 . 

Application of this method to image processing 
techniques demands a lot of computational effort. 
The present invention provides an implementation 
of an image processing method in which each pixel 
is replaced by a digital value that is obtained as a 
function of the original pixel value, the average 
value of pixel values within a window surrounding 
the pixel of interest and of the standard deviation of 
pixel values within such region. 
An exact value of the standard deviation of the 
pixel values within such window is obtained within 
a reasonable amount of computation time by the 
following elements that are characteristic of the 
present invention: 

-S? is obtained as the result of a convolution of the 
square values of the pixel values x of the image 
with a kernel having a number equal to N 2 .M 2 of 
identical elements surrounding pixel x that is pref- 
erably a power of 2, 

-said square values are obtained through the use of 
a look-up table, 

-the result of said convolution is divided by the 
number of pixels in said kernel namely N2.M2, 
-Si is obtained as the result of a second convolu- 
tion of the pixel values x of the image and a kernel 
having a number N1.M1 of identical elements sur- 
rounding pixel x that is preferably a power of 2, 
-the result of said convolution is divided by the 
number of pixels in said kernel, namely by Ni.Mi, 
-the square value (Si/(Ni .M1)) 2 of the result of the 
above-mentioned division is obtained through the 
use of a look-up table, 

-the subtraction S 2 /(N 2 .M 2 ) - (Sn/(Ni.Mi)) 2 is per- 
formed, 

-the square root of the result of the above subtrac- 
tion is obtained through the use of a look-up table. 

The method of the present invention comprises 
two convolutions performed on the pixel values or 
the square pixel values and kernels preferably 
comprising numbers of elements that are a power 
of 2. 

By using kernels having a number of elements 
that is a power of 2 the divisions needed to obtain 



an exact value of the standard deviation are re- 
duced to divisions by numbers that are a power of 
2. Normally divisions demand a lot of computa- 
tional effort. Now, the result of such divisions can 
5 be obtained through bit shifting which method de- 
mands less computational effort than methods re- 
quired to perform divisions by numbers that are not 
power of 2. 

The kernels may be square kernels so that Nt is 
io equal to Mi and N 2 is equal to M 2 . Preferably Ni, 
Mi , N 2 and M 2 are identical values. It is also 
possible that Nt is equal to N 2 and Mi is equal to 
M 2 . 

Further reduction of the computational effort is 

75 obtained by the frequent use of look-up tables. 

Once the look-up tables used to obtain square 
values and square root values are built and stored 
into the memory of an image processing computer, 
they are available for each image processing ap- 

20 piication.The look-up table for obtaining the square 
root values can e.g. be construed as follows: since 
for example sqrt(z . 2 10 ) = sqrt(z).2 5 , the calculation 
of the square root of a 20 bit number can be 
facilitated by using a look-up table for looking up 

25 the square root by means of the 10 most signifi- 

cant bits. If for example the look-up table is de- 

noted by R(z), z being a variable then R(z) is equal 
to sqrt(z).2 5 ; z = 0...1023. 

Alternatively S 2 can be obtained as the result 

30 of a convolution of a kernel of N 2 times M 2 iden- 
tical elements and values L(x) obtained through the 
use of a look-up table comprising the values of a 
substantially monotonous non-decreasing function. 
Preferably this function is an approximation of the 

35 function x 2 .AIso the second look-up table used to 
look-up the square root values may comprise a 
substantially monotonous non-decreasing function . 
Preferably this function is an approximation of the 
square root function. 

40 A particular aspect of the present invention 

relates to the adaptive histogram equalisation im- 
age processing technique. 

As has been described hereinbefore a prejudice 
exists against the incorporation of the standard 

45 deviation of the pixels in a local area around a pixel 
of interest because of the computation effort that is 
needed to exert such a method. 

it has been found through research and inves- 
tigation that with a certain illumination dose, giving 

so a particular average pixel value, the average histo- 
gram calculated by means of the histograms be- 
longing to local windows of an image (dimensions 
of a window being defined) can be approximated 
by a Gauss distribution with variable x, mean value 

55 m and standard deviation s (or standardized normal 
distribution with variable (x-m)/s). 
As has been explained hereinbefore, adaptive 
histogram equalisation is performed by computing 
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for each pixel a local histogram (in a pixels win- 
dow) and a rescaled cumulative histogram. Next, 
the pixel of interest is mapped to its corresponding 
value in the rescaled cumulative histogram. Taking 
into account the above research results that the 
average local histogram of the image areas cor- 
responding with a given value of the average pixel 
value m can be approximated by a normal distribu- 
tion, the cumulative histogram of which is known in 
the art, adaptive histogram equalisation can be i 
performed without the need of calculating the exact 
histograms of the image area. The resulting image 
after adaptive histogram equalisation processing 
taking into account the standard deviation can be 
formulated as follows: y = A . G((x-m)/s), ? 
wherein y is the resulting pixel value after image 
processing is performed, x is an original pixel val- 
, ue, A represents a scale factor, G represents the 
s standardized cumulative normal distribution, m the 
' average value and s is the standard deviation. 2c 
> Preferably m is obtained as St (Nt.Mi) wherein Si 
.Ni .and Mi are as defined hereinbefore. N,.Mi is 
preferably a power of 2. 

In an alternative embodiment any substantially mo- 
notonous non-decreasing function can be used in- 25 
stead of G . The scale factor A may be function of 
the pixel value x and/or the average value m and/or 
the standard deviation s. 

In a preferred embodiment of adaptive histo- 
gram equalisation , selected because of the mini- 30 
mum of computation time that is needed for cal- 
culation of the resulting image after image process- 
ing, the standard deviation is obtained according to 
the following formula 

s = (S 2 .'(N 2 .M2) - (Si/(Ni.M,)) 2 ) 1/2 as hereinabove 35 
defined wherein each of Si ,Ni ,Mi .S2.N2.and M 2 
have any of the definitions given hereinbefore with 
the preferences given and wherein the square root 
values are obtained through the use of a look-up 
table containing an approximation of the square 40 
root values. 

Further decrease of the computational effort is 
obtained in a preferred embodiment of the present 
invention wherein a (two-dimensional) look-up table 
is used for looking up the result of the formula A . 45 
G((x-m)/s), the inputs to such a look-up table are 
the result of a subtraction of the pixel value x and 
the average value m (wherein said value m is 
obtained as Si/(Ni.Mi) as defined hereinbefore 
)and the standard deviation s. 50 

Still further decrease of the computational effort 
can be obtained by combining the square root 
look-up table and the above mentioned two dimen- 
sional look-up table. 

The standard deviation is preferably obtained 55 
as is hereinabove described. This allows to reduce 
the computational effort to 14 additions and 5 bit- 
shifts per pixel. 



Alternatively the look-up table may comprise 
the result of the following formula: y = A.G((x-m)- 
-max(s. s min )) wherein s mm is a predetermined mini- 
mum value. This formula provides an alternative for 
5 ' the contrast limited adaptive histogram equalisation 
technique. 

G represents the standardized cumulative normal 
distribution. 

Alternatively G may be replaced by any monoto- 
0 nous non-decreasing function. 

Another particular aspect of the present inven- 
tion relates to the unsharp masking image process- 
ing technique. 

It was found that a very good result after image 
5 processing can be obtained by performing the fol- 
lowing transformation on the pixel values: 
y = x + B.g((x-m)/s) or-what is nearly equivalent - 
y = m + B.g((x-m)/s) ; 
wherein y represents a resulting pixel value after 
) image processing, x represents an original pixel 
value, B determines the amplification, m is an j 
average pixel v alue an d ^s^the-stano^rd^^ ? 
of tfie pixeTs^irTa region^sj^reundfng^^ / 
in^e^>cr^ie^function g(z) equals G(z) -1/2 / 
wherein ~G*(z) is the standardized cumulative Gaus- 
sian curve. Alternatively the function g.may be 
replaced by any function having the characteristic 
that g' is greater than or equal to -1 (g' represent- 
ing the derivate function). 

The scale factor B determining the amplification 
may be function of the pixel value x and/or the 
average value m and/or the standard deviation s. 
Commonly unsharp masking is followed by a 
smoothing step for the purpose of eliminating un- 
desirable effects (such as excessive noise en- 
hancement) resulting from said unsharp masking 
image processing. When performing the amended 
version of the unsharp masking technique accord- 
ing to the present invention wherein y = x + B.g(- 
(x-m)/s) or y = m + B.g((x-m)/s) said smoothing 
step may be omitted. 

In a preferred embodiment of the unsharp 
masking image processing , the standard deviation 
is obtained according to the formula described 
above . 

Decrease of the computational effort is ob- 
tained in a preferred embodiment of the present 
invention wherein a two-dimensional look-up table 
is used for looking up the result of the formula B . 
g((x-m)/s).The inputs to such a look-up table are 
the result of a subtraction of the pixel value x and 
the average value m (wherein said value m is 
obtained as Si/(N,.Mi) as defined above ) and the 
standard deviation s.The standard deviation is pref- 
erably obtained as is hereinabove described. 

Alternatively the two-dimensional look-up table 
may comprise the result of the following formula: B 
. g((x-m)/max(s, s mi n)) wherein s min is a predeter- 
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mined minimum value. 

As has already been explained hereinbefore, 
further decrease of the computational effort can be 
obtained by combining the used look-up tables. 

Claims 

1. An image processing method wherein a re- 
presentation of an image in the form of digital pixel 
values x is processed to obtain resulting digital 
pixel values y as a function of a pixel value x, an 
average value m of pixel values in a window sur- 
rounding said pixel x and a standard deviation s of 
pixel values in such window, characterised in that 
said standard deviation s is obtained as 
<S 2 /(N 2 .M 2 ) - (Si;(Ni.M,)) 2 ) ,/2 wherein S, is ob- 
tained as the result of a convolution of pixel values 
x and a kernel of Ni times Mi identical elements 
surrounding said pixel x and S 2 is the result of a 
convolution of a kernel of N 2 times M2 identical 
elements surrounding said pixel x and values L(x) 
obtained through the use of a first look-up table 
comprising the squares or the values of a substan- 
tially monotonous non-decreasing function and 
- wherein- - the - square" -root" valaes~are""obTatned " 
through the use of a second look-up table contain- 
ing an approximation of the square root values . 

2. An image processing method according to 
claim 1 wherein the values N1.M1 and/or N 2 .M 2 are 
a power of 2. 

3. An image processing method according to 
claim 1 or claim 2 wherein said second look-up 
table represents a substantially monotonous non- 
decreasing function. 

4. An image processing method wherein a re- 
presentation of an image in the form of digital pixel 
values x is processed to obtain resulting digital 
pixel values y as a function of a pixel value x, an 
average value m of pixel values in a window sur- 
rounding said pixel x, and a standard deviation s of 
pixel values in such window, characterised in that 
said resulting pixels y are obtained as A.G((x-m)- 
/max(s, Smin)) wherein A is a scale factor, G repre- 
sents a substantially monotonous non-decreasing 
function and s min is a predetermined minimum val- 
ue. 

5. An image processing method according to 
claim 4 wherein s min is 0 so that said resulting 
pixels y are obtained as A . G((x-m)/s) wherein A is 
a scale factor.and G represents a substantially mo- 
notonous non-decreasing function. 

6. A method according to claim 4 or 5 wherein 
G represents the standardized cumulative normal 
distribution. 

7. An image processing method according to 
any of claims 4 to 6 wherein A is a function of x 
and/or m and/or s. 



8. A method according to any of claims 4 to 7 
wherein the result of A.G((x-m)/max(s, s min )) is ob- 
tained through the use of a look-up table. 

9. An image processing method wherein a re- 
5 presentation of an image in the form of digital pixel 

values x is processed to obtain resulting digital 
pixel values y as a function of a pixel value x, an 
average value m of pixel values in a window sur- 
rounding said pixel x, and a standard deviation s of 

10 pixel values in such window, characterised in that 
said resulting pixels y are obtained as x + B . g((x- 
m)/max(s, s min )) or as m + B . g((x-m)/max(s, s mm )) 
wherein B determines the amplification ,s min is a 
predetermined minimum value and g(z) represents 

75 a function having a derivate that is equal to or 
larger than -1. 

10. An image processing method according to 
claim 9 wherein s min is 0 so that said resulting 
pixels y are obtained as x + B . g((x-m)/s) or as m 

20 + B . g((x-m)/s) wherein B determines the am- 
plification and g(z) represents a function having a 
derivate that is equal to or larger than -1 . 

11. A method according to claim 9 or 10 
wherein g(z) is a function equal to G(z) - 1/2 

25 wherein G is the standardized cumulative normal 
distribution: 

12. A method according to any of claims 9 to 
1 1 wherein B is a function of x and/or m and/or s. 

13. A method according to any of claims 9 to 
30 12 wherein the result of B .gUx-mymaxfs.Smin)) is 

obtained through the use of a two-dimentional look- 
up table. 

14. A method according to any of claims 4 to 
13 wherein s is obtained as ( S 2 /(N 2 .M 2 ) - 

35 (Si/(Ni.Mi)) 2 ) 1/2 wherein Si is obtained as the 
result of a convolution of pixel values x and a 
kernel of N1 times Mi identical elements surround- 
ing pixel value x and S 2 is the result of a convolu- 
tion of a kernel of N 2 times M 2 identical elements 

40 surrounding said pixel x and values L(x) obtained 
through the use of a first look-up table comprising 
the squares or the values of a substantially monot- 
onous non-decreasing function and wherein the 
square root values are obtained through the use of 

45 a second look-up table containing an approximation 
of the square root values. 

15. An image processing method according to 
claim 14 wherein the values Ni.Mi and/or N 2 .M 2 
are a power of 2. 

so 16. An image processing method according to 

claim 14 or 15 wherein said second look-up table 
represents a substantially monotonous non-de- 
creasing function. 

17. A method according to any of the preced- 

55 ing claims wherein said average value m is ob- 
tained as Si /N1.M1 wherein Si is the result of a 
convolution of the original image x with a kernel of 
N1 times Mi identical elements surrounding pixel 
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